{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d405c15d-533e-4a3c-a9a9-25792f41a633",
   "metadata": {},
   "source": [
    "### Enrichment of misexpression-associated SVs in regulatory annotations \n",
    "\n",
    "**Features**\n",
    "* 3D genome architecture \n",
    "    * TAD boundaries from: \n",
    "        * GM12878 (shared across multiple cell types) \n",
    "    * A/B compartments \n",
    "        * GM12878 cell line \n",
    "* CTCF-binding sites \n",
    "    * cCREs CTCF: all cell-types \n",
    "        - CTCF-only elements from all cell-types \n",
    "    * Specific primary cells: \n",
    "        * CD14 monocytes \n",
    "        * B-cells \n",
    "        * Neutrophils \n",
    "* Chromatin features \n",
    "    * Chrom HMM states \n",
    "* CpG islands \n",
    "\n",
    "--- \n",
    "\n",
    "* Logistic regression SV length adjusted and non SV length adjusted. \n",
    "* Model: misexpression-associated(1)/control(0) ~ feature (+ SV length)\n",
    "* Test duplications and deletions separately "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b70790a6-1632-48e0-a9db-2207a620ddec",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd \n",
    "import numpy as np\n",
    "from pybedtools import BedTool\n",
    "from io import StringIO\n",
    "from statsmodels.discrete import discrete_model\n",
    "from patsy import dmatrices\n",
    "from pathlib import Path\n",
    "from functools import reduce\n",
    "from statsmodels.stats import multitest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "44ee9daa",
   "metadata": {},
   "outputs": [],
   "source": [
    "wkdir = \"/lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/misexpression_v3/\"\n",
    "wkdir_path = Path(wkdir)\n",
    "# inputs \n",
    "sv_info_path = \"/lustre/scratch126/humgen/projects/interval_rna/interval_rna_seq/thomasVDS/lof_missense/data/sv_vcf/info_table/final_sites_critical_info_allele.txt\"\n",
    "all_vrnts_path = wkdir_path.joinpath(\"5_misexp_vrnts/test_cntrl_sets/vrnt_id_in_window_cntrl_misexp_genes.txt\")\n",
    "all_vrnts_bed_path = wkdir_path.joinpath(\"5_misexp_vrnts/test_cntrl_sets/vrnt_id_in_windows_misexp_genes.bed\")\n",
    "misexp_vrnts_path = wkdir_path.joinpath(\"5_misexp_vrnts/test_cntrl_sets/vrnt_id_misexp_tpm_zscore_median.txt\")\n",
    "# shared TAD boundaries across cell lines \n",
    "shared_tad_boundaries_path = wkdir_path.joinpath(\"reference/4d_nucleome/shared_boundaries/4DNFIVK5JOFU_imr90_huvec_hnek_hmec.bed\")\n",
    "# ENCODE cCREs \n",
    "encode_c_cres_dir = wkdir_path.joinpath(\"reference/encode/encode_c_cres\")\n",
    "# CpG islands \n",
    "cpg_isl_bed_path = wkdir_path.joinpath(\"reference/cpg_islands/cpgIslandExt.bed\")\n",
    "# PBMCs chromHMM \n",
    "pbmcs_chromhmm_15marks_path = wkdir_path.joinpath(\"reference/chromhmm/E062_15_coreMarks_hg38lift_mnemonics.bed\")\n",
    "# A/B compartments in GM12878\n",
    "gm12878_ab_path = wkdir_path.joinpath(\"reference/4d_nucleome/gm12878_hi_c/compartments/4DNFILYQ1PAY.bg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9d195f98",
   "metadata": {},
   "outputs": [],
   "source": [
    "# output\n",
    "out_dir = wkdir_path.joinpath(\"5_misexp_vrnts/functional\")\n",
    "out_dir_path = Path(out_dir)\n",
    "out_dir_path.mkdir(parents=True, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "75eeff15",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of variant IDs in windows: 20255\n",
      "Number of unique misexpression-associated SVs: 105\n"
     ]
    }
   ],
   "source": [
    "# load bed file with variants in windows \n",
    "vrnts_in_windows_bed = BedTool(all_vrnts_bed_path)\n",
    "vrnts_in_windows_bed_sorted = vrnts_in_windows_bed.sort()\n",
    "# all variants in windows df\n",
    "all_vrnts_in_windows_df = pd.read_csv(all_vrnts_path, sep=\"\\t\", header=None).rename(columns={0:\"vrnt_id\"})\n",
    "all_vrnts_in_windows = all_vrnts_in_windows_df.vrnt_id.unique()\n",
    "print(f\"Number of variant IDs in windows: {len(all_vrnts_in_windows)}\")\n",
    "\n",
    "# load misexpression-associated variants \n",
    "misexp_vrnts_ids = pd.read_csv(misexp_vrnts_path, sep=\"\\t\", header=None)[0].astype(str).unique()\n",
    "print(f\"Number of unique misexpression-associated SVs: {len(misexp_vrnts_ids)}\")\n",
    "all_vrnts_in_windows_df[\"misexp_uniq\"] = np.where(all_vrnts_in_windows_df.vrnt_id.isin(misexp_vrnts_ids), 1, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "beef430e",
   "metadata": {},
   "source": [
    "### Shared TAD boundaries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1b4fad00-ad84-4098-a18d-d730d438c58b",
   "metadata": {},
   "outputs": [],
   "source": [
    "tad_bed_intersect_cols = {0:\"vrnt_chrom\", 1:\"vrnt_start\", 2:\"vrnt_end\", 3:\"vrnt_id\", 4:\"boundary_chrom\", \n",
    "                      5:\"boundary_start\", 6:\"boundary_end\", 7:\"boundary_no\", 8:\"boundary_score\", 9:\"overlap\"}\n",
    " \n",
    "vrnt_tad_features_df = all_vrnts_in_windows_df.copy()\n",
    "# TAD boundaries bed file \n",
    "tad_boundaries_path = shared_tad_boundaries_path\n",
    "tad_boundaries_bed = BedTool(tad_boundaries_path)\n",
    "input_tad_bed = tad_boundaries_bed.sort()\n",
    "\n",
    "# identify SVs that overlap TAD boundaries\n",
    "sv_intersect_tad_str = StringIO(str(vrnts_in_windows_bed_sorted.intersect(input_tad_bed, wo=True)))\n",
    "sv_intersect_tad_df = pd.read_csv(sv_intersect_tad_str, sep=\"\\t\", header=None).rename(columns=tad_bed_intersect_cols).astype({\"vrnt_id\":str})\n",
    "sv_intersect_tad = sv_intersect_tad_df.vrnt_id.unique()\n",
    "\n",
    "# annotate variants intersecting TAD boundaries and window around TAD\n",
    "vrnt_tad_features_df[f\"gm12878_shared_intersect_tad_boundary\"] = np.where(vrnt_tad_features_df.vrnt_id.isin(sv_intersect_tad), 1, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a32abfd5",
   "metadata": {},
   "source": [
    "### A/B compartments "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "017810fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load A/B compartments in GM12878\n",
    "gm12878_ab_bed = BedTool(gm12878_ab_path)\n",
    "gm12878_ab_input = gm12878_ab_bed.sort()\n",
    "\n",
    "# intersect with A/B compartment scores \n",
    "ab_bed_intersect_cols = {0: 'vrnt_chrom', 1: 'vrnt_start', 2: 'vrnt_end', 3: 'vrnt_id', \n",
    "                           4: 'compartment_chrom', 5: 'compartment_start',6: 'compartment_end',\n",
    "                           7: 'score', 8: 'overlap'}\n",
    "sv_intersect_ab_str = StringIO(str(vrnts_in_windows_bed_sorted.intersect(gm12878_ab_input, wo=True)))\n",
    "sv_intersect_ab_df = pd.read_csv(sv_intersect_ab_str, sep=\"\\t\", header=None).rename(columns=ab_bed_intersect_cols).astype({\"vrnt_id\":str})\n",
    "# label compartment types \n",
    "conditions = [(sv_intersect_ab_df.score >= 0) & (~sv_intersect_ab_df.score.isnull()), \n",
    "              (sv_intersect_ab_df.score < 0) & (~sv_intersect_ab_df.score.isnull()), \n",
    "              (sv_intersect_ab_df.score.isnull())]\n",
    "values = [\"A\", \"B\", \"Unassigned\"]\n",
    "\n",
    "sv_intersect_ab_df[\"compartment_type\"] = np.select(conditions, values)\n",
    "vrnt_compartment_features_df = all_vrnts_in_windows_df.copy()\n",
    "# annotate SV with binary variable with compartment overlap\n",
    "for overlap in [\"A\", \"B\"]:\n",
    "    vrnt_overlap_compartment = sv_intersect_ab_df[sv_intersect_ab_df.compartment_type == overlap].vrnt_id.unique()\n",
    "    vrnt_compartment_features_df[f\"{overlap}_overlap\"] = np.where(vrnt_compartment_features_df.vrnt_id.isin(vrnt_overlap_compartment), 1, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b790e13-bf60-4597-a21c-a4a637b11e70",
   "metadata": {},
   "source": [
    "### CTCF-bound cCREs "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "0dd138ed-2461-4efb-862b-2d7b60339d28",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "all: number of unique variants in dataframe: 20255\n",
      "CD14_monocyte: number of unique variants in dataframe: 20255\n",
      "B_cell: number of unique variants in dataframe: 20255\n",
      "Neutrophil: number of unique variants in dataframe: 20255\n"
     ]
    }
   ],
   "source": [
    "bed_intersect_cols = {0:\"vrnt_chrom\", 1:\"vrnt_start\", 2:\"vrnt_end\", 3:\"vrnt_id\"}\n",
    "encode_c_cres_dir_path = Path(encode_c_cres_dir)\n",
    "\n",
    "bed_path_list = [\"GRCh38-cCREs.CTCF-only.bed\", \n",
    "                 \"ENCFF389PZY_ENCFF587XGD_ENCFF184NWF_ENCFF496PSJ.7group.bed\", \n",
    "                 \"ENCFF035DJL.7group.bed\", \n",
    "                 \"ENCFF685DZI_ENCFF311TAY_ENCFF300LXQ.7group.bed\"\n",
    "                ]\n",
    "\n",
    "cell_types_list = [\"all\", \n",
    "                   \"CD14_monocyte\", \n",
    "                   \"B_cell\", \n",
    "                   \"Neutrophil\"\n",
    "                  ]\n",
    "\n",
    "vrnt_ccre_features_df = all_vrnts_in_windows_df.copy() \n",
    "for cell_type, bed_file in zip(cell_types_list, bed_path_list): \n",
    "    ctcf_peaks_path = encode_c_cres_dir_path.joinpath(f\"{bed_file}\")\n",
    "    peak_bed = BedTool(ctcf_peaks_path).sort()\n",
    "    all_chrom_peak_info = []   \n",
    "    # identify SVs that overlap peaks\n",
    "    sv_intersect_str = StringIO(str(vrnts_in_windows_bed_sorted.intersect(peak_bed, wo=True)))\n",
    "    # check if string is empty \n",
    "    if not sv_intersect_str.getvalue():\n",
    "        raise ValueError(\"No SVs overlap ChIP-seq peaks\")\n",
    "    sv_intersect_df = pd.read_csv(sv_intersect_str, sep=\"\\t\", header=None).astype({3:str}).rename(columns=bed_intersect_cols)\n",
    "    if cell_type == \"all\": \n",
    "        cre_type_index = -2\n",
    "    else: \n",
    "        cre_type_index = -3\n",
    "    cre_types_list = sv_intersect_df[sv_intersect_df.columns[cre_type_index]].unique().tolist()\n",
    "    for cre_type in cre_types_list: \n",
    "        vrnt_id_intersect_cre_type = sv_intersect_df[sv_intersect_df[sv_intersect_df.columns[cre_type_index]] == cre_type].copy().vrnt_id.unique()\n",
    "        vrnt_ccre_features_df[f\"{cre_type.replace('-', '').replace(',', '')}_{cell_type}\"] = np.where(vrnt_ccre_features_df.vrnt_id.isin(vrnt_id_intersect_cre_type), 1, 0)\n",
    "    print(f\"{cell_type}: number of unique variants in dataframe: {len(vrnt_ccre_features_df.vrnt_id.unique())}\")\n",
    "# subset to CTCF features \n",
    "vrnt_ctcf_ccre_features_df = vrnt_ccre_features_df[[\"vrnt_id\", \"misexp_uniq\", \"CTCFonlyCTCFbound_all\", 'CTCFonlyCTCFbound_CD14_monocyte', 'HighCTCF_B_cell', \"HighCTCF_Neutrophil\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba7ea739-b942-4f0d-b671-571ab13519cc",
   "metadata": {},
   "source": [
    "### CpG Islands "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b00bac3e-9d00-4651-9757-3a852a599780",
   "metadata": {},
   "outputs": [],
   "source": [
    "### CpG Islands \n",
    "cpg_isl_bed = BedTool(cpg_isl_bed_path).sort()\n",
    "# identify SVs that overlap islands\n",
    "bed_intersect_cols_cpg = {0:\"vrnt_chrom\", 1:\"vrnt_start\", 2:\"vrnt_end\", 3:\"vrnt_id\", 4:\"cpg_isl_chrom\", \n",
    "                          5:\"cpg_isl_start\", 6:\"cpg_isl_end\", 7:\"overlap\"}\n",
    "sv_intersect_cpg_isl_str = StringIO(str(vrnts_in_windows_bed_sorted.intersect(cpg_isl_bed, wo=True)))\n",
    "sv_intersect_cpg_isl_df = pd.read_csv(sv_intersect_cpg_isl_str, sep=\"\\t\", header=None).rename(columns=bed_intersect_cols_cpg).astype({\"vrnt_id\":str})\n",
    "sv_intersect_cpg_isl = sv_intersect_cpg_isl_df.vrnt_id.unique()\n",
    "# annotate variants intersecting at least one CpG island\n",
    "vrnts_cpg_island_df = all_vrnts_in_windows_df.copy() \n",
    "vrnts_cpg_island_df[\"intersect_cpg_isl\"] = np.where(vrnts_cpg_island_df.vrnt_id.isin(sv_intersect_cpg_isl), 1, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12ed990b-a22d-4dd2-8b85-30dc3b8edbf4",
   "metadata": {},
   "source": [
    "### Chromatin States "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "94146b1d-a5d5-4990-8b4e-037024d705e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "chromhmm_15states = ['1_TssA', '2_TssAFlnk', '3_TxFlnk', '4_Tx', '5_TxWk', '6_EnhG', '7_Enh', '8_ZNF/Rpts', '9_Het', \n",
    "                     '10_TssBiv', '11_BivFlnk', '12_EnhBiv', '13_ReprPC', '14_ReprPCWk', '15_Quies']\n",
    "\n",
    "bed_intersect_cols_chrom_hmm = {0:\"vrnt_chrom\", 1:\"vrnt_start\", 2:\"vrnt_end\", 3:\"vrnt_id\", 4:\"state_chrom\", \n",
    "                                5:\"state_start\", 6:\"state_end\", 7:\"state_name\", 8:\"overlap\"}\n",
    "\n",
    "# load PBMCs chrom HMM states\n",
    "pbmcs_chromhmm_15marks_bed = BedTool(pbmcs_chromhmm_15marks_path)\n",
    "    \n",
    "# identify SVs that overlap peaks\n",
    "sv_intersect_str = StringIO(str(vrnts_in_windows_bed_sorted.intersect(pbmcs_chromhmm_15marks_bed, wo=True)))\n",
    "sv_intersect_df = pd.read_csv(sv_intersect_str, sep=\"\\t\", header=None).astype({3:str}).rename(columns=bed_intersect_cols_chrom_hmm)\n",
    "# build dataframe \n",
    "vrnts_chrom_hmm_df = all_vrnts_in_windows_df.copy()\n",
    "for state in chromhmm_15states: \n",
    "    # check if variant overlaps state \n",
    "    vrnt_ids_with_state = sv_intersect_df[sv_intersect_df.state_name == state].vrnt_id\n",
    "    vrnts_chrom_hmm_df[state] = np.where(vrnts_chrom_hmm_df.vrnt_id.isin(set(vrnt_ids_with_state)), 1, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77236e1c-2094-4771-9e43-476c98ba1fd9",
   "metadata": {},
   "source": [
    "### Combine features "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "598f09c8-a6ff-44d8-a5df-fef4a84bd590",
   "metadata": {},
   "outputs": [],
   "source": [
    "# merge different features \n",
    "dfs_to_merge = [vrnt_tad_features_df, \n",
    "                vrnt_compartment_features_df, \n",
    "                vrnt_ctcf_ccre_features_df, \n",
    "                vrnts_cpg_island_df,\n",
    "                vrnts_chrom_hmm_df]\n",
    "\n",
    "vrnt_features_merged_df = reduce(lambda  left,right: pd.merge(left,right, on=all_vrnts_in_windows_df.columns.tolist(),\n",
    "                                                              how='inner'), dfs_to_merge)\n",
    "\n",
    "# remove / in name and remove number from name for chromatin HMM columns \n",
    "chromhmm_rename_columns = {col:\"\".join(col.split(\"/\")).split(\"_\")[1] for col in chromhmm_15states}\n",
    "vrnt_features_merged_df = vrnt_features_merged_df.rename(columns=chromhmm_rename_columns)\n",
    "if set(vrnt_features_merged_df.vrnt_id.unique()) != set(all_vrnts_in_windows): \n",
    "    raise ValueError(\"SVs with functional annotations does not match input variants.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d284a05f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# add structural variant information \n",
    "sv_info_df =pd.read_csv(sv_info_path, sep=\"\\t\", dtype={\"plinkID\": str}).rename(columns={\"plinkID\":\"vrnt_id\"})\n",
    "\n",
    "vrnt_features_merged_info_df = pd.merge(vrnt_features_merged_df, \n",
    "                                   sv_info_df, \n",
    "                                   on=\"vrnt_id\", \n",
    "                                   how=\"left\"\n",
    "                                  )\n",
    "# check for NaNs\n",
    "if vrnt_features_merged_info_df.isnull().values.any(): \n",
    "    raise ValueError(\"SV functional annotation dataframe contains NaNs.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "04a26b83",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create features directory\n",
    "features_dir_path = out_dir_path.joinpath(\"features\")\n",
    "features_dir_path.mkdir(parents=True, exist_ok=True)\n",
    "# write features to file \n",
    "vrnt_features_out = features_dir_path.joinpath(\"vrnt_features_reg_annot.csv\")\n",
    "vrnt_features_merged_info_df.to_csv(vrnt_features_out, index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a83fd1b3-5f29-4e5d-ba89-099184282f58",
   "metadata": {},
   "source": [
    "### Enrichment analysis \n",
    "\n",
    "* Normalise all features (including binary) and do single feature logistic regression "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "e5b32c64-aa09-4fac-8a65-15c692bf713c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# group features into 3D genome and epigenetic features \n",
    "features_dict = {\"3d_genome\": ['gm12878_shared_intersect_tad_boundary', 'A_overlap', 'B_overlap', 'CTCFonlyCTCFbound_all', \n",
    "                                'CTCFonlyCTCFbound_CD14_monocyte', 'HighCTCF_B_cell', \"HighCTCF_Neutrophil\",\n",
    "                              ], \n",
    "                 \"epigenetics\":[\"intersect_cpg_isl\"] + list(chromhmm_rename_columns.values())\n",
    "                }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "e3be963c-4927-40ee-a117-259dbf282acc",
   "metadata": {
    "jupyter": {
     "outputs_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030222\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030460\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030443\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030402\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030368\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030115\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030180\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030096\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030312\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030462\n",
      "         Iterations 9\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/nfs/users/nfs_t/tv5/.conda/envs/tv5_base/lib/python3.7/site-packages/statsmodels/discrete/discrete_model.py:1819: RuntimeWarning: overflow encountered in exp\n",
      "  return 1/(1+np.exp(-X))\n",
      "/nfs/users/nfs_t/tv5/.conda/envs/tv5_base/lib/python3.7/site-packages/statsmodels/discrete/discrete_model.py:1872: RuntimeWarning: divide by zero encountered in log\n",
      "  return np.sum(np.log(self.cdf(q*np.dot(X,params))))\n",
      "/nfs/users/nfs_t/tv5/.conda/envs/tv5_base/lib/python3.7/site-packages/statsmodels/base/model.py:606: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  ConvergenceWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: Maximum number of iterations has been exceeded.\n",
      "         Current function value: inf\n",
      "         Iterations: 1000\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.029823\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030164\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030375\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030172\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030454\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030413\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030450\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030378\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030092\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030327\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030166\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030425\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.049431\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.047650\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.048130\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.049553\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.048034\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.048253\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.047773\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.046654\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.044903\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.045504\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.048782\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.046832\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.043372\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.045932\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.042952\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.050338\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.049756\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.049472\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.050018\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.047638\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.046920\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.048244\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.049936\n",
      "         Iterations 10\n"
     ]
    }
   ],
   "source": [
    "# logistic regression \n",
    "logistic_regr_results, logr_count = {}, 0         \n",
    "for sv_type in [\"DEL\", \"DUP\"]:\n",
    "    input_df = vrnt_features_merged_info_df[vrnt_features_merged_info_df.SVTYPE == sv_type].copy()\n",
    "    input_df[\"svlen_norm\"] = (input_df[\"SVLEN\"] - input_df[\"SVLEN\"].mean())/input_df[\"SVLEN\"].std()\n",
    "    for key in features_dict.keys():\n",
    "        features = features_dict[key]\n",
    "        for feature in features:\n",
    "            input_df[f\"{feature}_norm\"] = (input_df[feature] - input_df[feature].mean())/input_df[feature].std()\n",
    "            y, X = dmatrices(f'misexp_uniq ~ {feature}_norm + svlen_norm', input_df, return_type = 'dataframe')\n",
    "            logit_fit = discrete_model.Logit(endog=y, exog=X).fit(maxiter=1000)\n",
    "            log_odds, pval = logit_fit.params[1], logit_fit.pvalues[1]\n",
    "            # normal approximation confidence intervals\n",
    "            lower_conf = logit_fit.conf_int(alpha=0.05)[0][1]\n",
    "            upper_conf = logit_fit.conf_int(alpha=0.05)[1][1]\n",
    "            logistic_regr_results[logr_count] = [key, feature, sv_type, log_odds, lower_conf, upper_conf, pval]\n",
    "            logr_count += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "1790e119-70d9-48f6-9eff-43e639cc3d56",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "columns_logr_enrich = [\"category\", \"feature\", \"sv_type\", \"log_odds\", \"lower\", \"upper\", \"pval\"]\n",
    "logistic_regr_results_df = pd.DataFrame.from_dict(logistic_regr_results, orient=\"index\", columns=columns_logr_enrich)\n",
    "\n",
    "# multiple testing correction\n",
    "pval_list = logistic_regr_results_df.pval.to_numpy()\n",
    "# BH FDR\n",
    "reject, pvals_corrected, _, _ = multitest.multipletests(pval_list, alpha=0.05, method=\"fdr_bh\")\n",
    "logistic_regr_results_df[\"pass_fdr_bh\"] = reject\n",
    "logistic_regr_results_df[\"pvals_corrected_fdr_bh\"] = pvals_corrected\n",
    "# Bonferroni correction \n",
    "reject, pvals_corrected, _, _ = multitest.multipletests(pval_list, alpha=0.05, method=\"bonferroni\")\n",
    "logistic_regr_results_df[\"pass_bonf\"] = reject\n",
    "logistic_regr_results_df[\"pvals_corrected_bonf\"] = pvals_corrected"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a68bd17d",
   "metadata": {},
   "source": [
    "### Significant results "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "9fd89cbd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enh DUP 1.2620667347929995\n",
      "TxWk DUP 1.2483521290660107\n",
      "intersect_cpg_isl DUP 0.8790587581418776\n",
      "TssA DUP 0.8279293576899488\n",
      "TssAFlnk DUP 0.8131586189287672\n",
      "ReprPC DUP 0.7166115443339292\n",
      "Tx DUP 0.6961995041587604\n",
      "EnhG DUP 0.5916922200004098\n",
      "Tx DEL 0.36268782325369375\n",
      "ReprPCWk DEL 0.3079152534175603\n",
      "HighCTCF_B_cell DEL 0.297143856750386\n",
      "TxWk DEL 0.2958800020260846\n",
      "intersect_cpg_isl DEL 0.27403255506016205\n",
      "Enh DEL 0.2504021197975904\n",
      "HighCTCF_Neutrophil DEL 0.23921334073247177\n",
      "gm12878_shared_intersect_tad_boundary DEL 0.19924272289615447\n",
      "EnhBiv DEL 0.19562396283062802\n"
     ]
    }
   ],
   "source": [
    "logistic_regr_results_pass_bonf_df = logistic_regr_results_df[logistic_regr_results_df.pass_bonf].sort_values(by=\"log_odds\", ascending=False)\n",
    "for index, row in logistic_regr_results_pass_bonf_df.iterrows(): \n",
    "    print(row[\"feature\"], row[\"sv_type\"], row[\"log_odds\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5420509",
   "metadata": {},
   "source": [
    "### Write results to file "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "ed966d09",
   "metadata": {},
   "outputs": [],
   "source": [
    "# category names \n",
    "category_names = {\"3d_genome\": \"3D Genome\", \"epigenetics\": \"Regulatory\"}\n",
    "# feature names \n",
    "feature_names = {'gm12878_shared_intersect_tad_boundary': \"TAD boundaries (shared)\", \n",
    "                 'CTCFonlyCTCFbound_all': \"CTCF (all)\",\n",
    "                 'CTCFonlyCTCFbound_CD14_monocyte': \"CTCF (CD14+ monocytes)\", \n",
    "                 'HighCTCF_B_cell': \"CTCF (B-cells)\",\n",
    "                 'HighCTCF_Neutrophil': \"CTCF (Neutrophils)\", \n",
    "                 'intersect_cpg_isl': \"CpG islands\", \n",
    "                 'TssA': \"Active TSS\", \n",
    "                 'TssAFlnk': \"Flanking active TSS\",\n",
    "                 'TxFlnk': \"Transcr. at gene 5' and 3'\", \n",
    "                 'Tx': \"Strong transcription\", \n",
    "                 'TxWk': \"Weak transcription\", \n",
    "                 'EnhG': \"Genic enhancers\", \n",
    "                 'Enh': \"Enhancers\", \n",
    "                 'ZNFRpts': \"ZNF genes & repeats\", \n",
    "                 'Het': \"Heterochromatin\", \n",
    "                 'TssBiv': \"Bivalent/poised TSS\",\n",
    "                 'BivFlnk': \"Flanking bivalent TSS/enhancer\", \n",
    "                 'EnhBiv': \"Bivalent enhancer\", \n",
    "                 'ReprPC': \"Repressed polyComb\", \n",
    "                 'ReprPCWk': \"Weak repressed polycomb\", \n",
    "                 'Quies': \"Quiescent\", \n",
    "                 'A_overlap': \"A compartment\", \n",
    "                 'B_overlap': \"B compartment\"\n",
    "                }\n",
    "logistic_regr_results_df[\"feature_name\"] = logistic_regr_results_df.feature.replace(feature_names)\n",
    "logistic_regr_results_df[\"category_name\"] = logistic_regr_results_df.category.replace(category_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "7c56dda8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# all results for deletions and duplications, use Bonferroni cutoff\n",
    "logistic_regr_results_df[\"log_odds_adj\"] = np.where(logistic_regr_results_df.pass_bonf, \n",
    "                                                    logistic_regr_results_df.log_odds, \n",
    "                                                    np.nan)\n",
    "# create results directory\n",
    "results_dir_path = out_dir_path.joinpath(\"results\")\n",
    "results_dir_path.mkdir(parents=True, exist_ok=True)\n",
    "# write results to file\n",
    "all_results_path = results_dir_path.joinpath(\"logr_func_enrich_results_all.tsv\")\n",
    "logistic_regr_results_df.to_csv(all_results_path, sep=\"\\t\", index=False)                                      "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20e84e06",
   "metadata": {},
   "source": [
    "### Non-length adjusted enrichments "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "d7e4be71",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030299\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030644\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030633\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030489\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030485\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030172\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030239\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030176\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030417\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030625\n",
      "         Iterations 9\n",
      "Warning: Maximum number of iterations has been exceeded.\n",
      "         Current function value: 0.030639\n",
      "         Iterations: 1000\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.029946\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030316\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030505\n",
      "         Iterations 9\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/nfs/users/nfs_t/tv5/.conda/envs/tv5_base/lib/python3.7/site-packages/statsmodels/base/model.py:606: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  ConvergenceWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030278\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030609\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030513\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030583\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030490\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030164\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030430\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030296\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.030624\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.050388\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.050647\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.051551\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.050535\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.048839\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.048971\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.048330\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.047158\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.045516\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.046238\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.051113\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.048057\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.044578\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.046893\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.043429\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.052576\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.050928\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.050554\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.051687\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.048430\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.047370\n",
      "         Iterations 9\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.049510\n",
      "         Iterations 10\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.052227\n",
      "         Iterations 10\n"
     ]
    }
   ],
   "source": [
    "# logistic regression \n",
    "logistic_regr_results_non_len_adj, logr_count = {}, 0\n",
    "for sv_type in [\"DEL\", \"DUP\"]:\n",
    "    input_df = vrnt_features_merged_info_df[vrnt_features_merged_info_df.SVTYPE == sv_type].copy()\n",
    "    for key in features_dict.keys():\n",
    "        features = features_dict[key]\n",
    "        for feature in features:\n",
    "            input_df[f\"{feature}_norm\"] = (input_df[feature] - input_df[feature].mean())/input_df[feature].std()\n",
    "            y, X = dmatrices(f'misexp_uniq ~ {feature}_norm', input_df, return_type = 'dataframe')\n",
    "            logit_fit = discrete_model.Logit(endog=y, exog=X).fit(maxiter=1000)\n",
    "            log_odds, pval = logit_fit.params[1], logit_fit.pvalues[1]\n",
    "            # normal approximation confidence intervals\n",
    "            lower_conf = logit_fit.conf_int(alpha=0.05)[0][1]\n",
    "            upper_conf = logit_fit.conf_int(alpha=0.05)[1][1]\n",
    "            logistic_regr_results_non_len_adj[logr_count] = [key, feature, sv_type, log_odds, lower_conf, upper_conf, pval]\n",
    "            logr_count += 1\n",
    "            \n",
    "logistic_regr_results_non_len_adj_df = pd.DataFrame.from_dict(logistic_regr_results_non_len_adj, orient=\"index\", columns=columns_logr_enrich)\n",
    "\n",
    "# multiple testing correction\n",
    "pval_list = logistic_regr_results_non_len_adj_df.pval.to_numpy()\n",
    "# BH FDR\n",
    "reject, pvals_corrected, _, _ = multitest.multipletests(pval_list, alpha=0.05, method=\"fdr_bh\")\n",
    "logistic_regr_results_non_len_adj_df[\"pass_fdr_bh\"] = reject\n",
    "logistic_regr_results_non_len_adj_df[\"pvals_corrected_fdr_bh\"] = pvals_corrected\n",
    "# Bonferroni correction \n",
    "reject, pvals_corrected, _, _ = multitest.multipletests(pval_list, alpha=0.05, method=\"bonferroni\")\n",
    "logistic_regr_results_non_len_adj_df[\"pass_bonf\"] = reject\n",
    "logistic_regr_results_non_len_adj_df[\"pvals_corrected_bonf\"] = pvals_corrected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "de98f9d3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enh DUP 1.3483934099463182\n",
      "TxWk DUP 1.3371875314903867\n",
      "intersect_cpg_isl DUP 0.9900244119091235\n",
      "TssA DUP 0.9199302894668132\n",
      "TssAFlnk DUP 0.9146454777560912\n",
      "HighCTCF_B_cell DUP 0.9067217923366605\n",
      "ReprPC DUP 0.8180668333601342\n",
      "Tx DUP 0.8038045067063008\n",
      "HighCTCF_Neutrophil DUP 0.7961936907052467\n",
      "CTCFonlyCTCFbound_CD14_monocyte DUP 0.7197906823358683\n",
      "EnhG DUP 0.6773341661297569\n",
      "EnhBiv DUP 0.6090085756204545\n",
      "gm12878_shared_intersect_tad_boundary DUP 0.5217493210617609\n",
      "TssBiv DUP 0.47056676524008567\n",
      "Tx DEL 0.37712889316068826\n",
      "HighCTCF_B_cell DEL 0.3295353732262771\n",
      "ReprPCWk DEL 0.3292197709590169\n",
      "TxWk DEL 0.30888578279640033\n",
      "intersect_cpg_isl DEL 0.3007085847169951\n",
      "Enh DEL 0.2767984206547811\n",
      "TxFlnk DUP 0.27517305164021144\n",
      "HighCTCF_Neutrophil DEL 0.2701904414729294\n",
      "gm12878_shared_intersect_tad_boundary DEL 0.22690569149468284\n",
      "EnhBiv DEL 0.21838755116601158\n",
      "TssA DEL 0.20067859853742795\n"
     ]
    }
   ],
   "source": [
    "logistic_regr_results_non_len_adj_pass_bonf_df = logistic_regr_results_non_len_adj_df[logistic_regr_results_non_len_adj_df.pass_bonf].sort_values(by=\"log_odds\", ascending=False)\n",
    "for index, row in logistic_regr_results_non_len_adj_pass_bonf_df.iterrows(): \n",
    "    print(row[\"feature\"], row[\"sv_type\"], row[\"log_odds\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "a12a8e82",
   "metadata": {},
   "outputs": [],
   "source": [
    "# add feature and category names\n",
    "logistic_regr_results_non_len_adj_df[\"feature_name\"] = logistic_regr_results_non_len_adj_df.feature.replace(feature_names)\n",
    "logistic_regr_results_non_len_adj_df[\"category_name\"] = logistic_regr_results_non_len_adj_df.category.replace(category_names)\n",
    "\n",
    "# all results for deletions and duplications, use Bonferroni cutoff\n",
    "logistic_regr_results_non_len_adj_df[\"log_odds_adj\"] = np.where(logistic_regr_results_non_len_adj_df.pass_bonf, \n",
    "                                                                logistic_regr_results_non_len_adj_df.log_odds, \n",
    "                                                                np.nan)\n",
    "all_results_no_len_adj_path = results_dir_path.joinpath(\"logr_func_enrich_results_all_no_len_adj.tsv\")\n",
    "logistic_regr_results_non_len_adj_df.to_csv(all_results_no_len_adj_path, sep=\"\\t\", index=False)   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9df80293",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
